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Abstract 

In this paper, the realization-free model approximation problem, as stated in 
dia, is revisited in the case where the interpolating model might he time-delay 
dependent. To this aim, the Loewner framework, initially settled for delay-free 
realization, is firstly generalized to the single delay case. Secondly, the (infinite) 
model approximation ,^2 optimality conditions are established through the use of 
the Lambert functions. Finally, a numerically effective iterative scheme, named 
dTF-IRKA, similar to the TF-IRKA O, is proposed to reach a part of the afore¬ 
mentioned optimality conditions. The proposed method validity and interest are 
assessed on different numerical examples. 


1 Introduction 

1.1 Motivating context and problem formulation 

Mathematical dynamical models are usually necessary to understand, analyse and con¬ 
trol the behaviour of physical phenomena. Usually, high fidelity models require numer¬ 
ous equations and variables. The resulting associated state-space realization is conse¬ 
quently of large dimensions and the resulting system is said to be a large-scale one. In 
addition, in some cases, a linear finite dimension realization is not accessible or even 
does not exist (e.g., partial differential equations models, irrational transfer functions, 
etc.). Although such models can faithfully and accurately reproduce reality, it might 
lead to (i) high complexity and/or (//j realization-less models for which classical meth¬ 
ods cannot be reasonably applied due to high numerical burden, low computational 
speed and inappropriate tools. In these cases and for control concerns, an approxima¬ 
tion by a less complex realization is therefore desirable. This justifies the use of model 
approximation and interpolation techniques allowing to find a simpler model which 
faithfully approaches the original one and that can be used in place for simulation, 
analysis and control (for survey and historical references, see Eiiiia and references 
therein). 
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Moreover, as Time-Delay Systems (TDS) is a large class of dynamical systems 
which generalizes the finite dimension realization one, approximating any transfer 
G or complex data set^(i,',i/(i,)) (for i = I,... ,r, r G N*), 

by a time-delay dynamical model might be relevant for some specific applications 
where the delay naturally appears. Indeed for such kind of systems many dedicated 
and powerful results have been obtained for stability, performance analysis and control 
(see e.g., EUlia). 

In this paper, the approximation of any realization or realization-free linear dynam¬ 
ical model by a single time-delay model of finite dimension is developed. More specif¬ 
ically, we are interested in approximating any MIMO transfer function H{s) G 
by a single delay finite-dimensional linear time-invariant descriptor system denoted 
Hd = (E,A,B,C,t) and defined by: 

£x(t) =Ax(t — t) -GBuf), y(t) = Cx{t), (1) 

whose transfer function is Hj{s) = C{sE — It is straightforward to note 

that the approximation form Q generalizes the delay-free one used in Ul |2l given as 
H= {E,A,B,C,0) (or simply {E,A,B,C)), 

Ex(t) =Ax(t) -\-Bu{t), y(t) = Cx(f). (2) 

Following lUJIll, and inspired by the widely used J^^-approximation problem ifTOl 
[mEi, our objective can be mathematically formulated as follows: 

Problem 1 Given a LTI system H{s) G (or {^Si,El{sifj, the evaluation of H(s) at 
Si G C, for i = l,...,r), a positive integer r G N* and a positive scalar T G K, find a 
model Hd = {E,A,B,C,t) G such that 

Hd:= argmin ||H-Gd||,^r 2 - (3) 

GdG.^.dim(Gd)<)' 


In other words, if an evaluation of the transfer function H{s), for any s G C, is 
available (either from data or by simply evaluating H{s)), our goal is to find a delay 
model of the form Q, that well approximates El, in the sense of the J^-norm. 

1.2 Contributions 

The purpose of this paper is thus to extend the application domain of the Loewner 
framework established in iniiia to dynamical systems with one single internal delay. 
To do so, a new delay Loewner framework is firstly developed to interpolate a given 
transfer function by a time-delay model of the form 0, enabling the delayed Loewner 
framework to be applied to any models for which the transfer function is accessible 
only. This allows then model approximation for both infinite/finite dimensional sys¬ 
tems and data-based ones. Then, following Problem the J^ 2 -oriented optimality 
conditions are formulated and used to construct an iterative algorithm, similar to the 
recently proposed TF-IRKA |j2), allowing to obtain an approximated model Hd satis¬ 
fying a finite number of the optimality conditions. 

'Throughout this paper, we denote or simply M 2 , the open subspace of M 2 with matrix-valued 

function H(s) (riy outputs, inputs), Vx G C, which are analytic in Re(i) > 0 (functions that are locally 
given by a convergent power series and differentiable on each point of its definition set) 0. 

^We denote as (^Si,H{si)) the evaluation of transfer H at Sj. 



1.3 Notations and outlines 


We denote by N* the set of natural numbers without 0, by M 2 the Hilbert space of 
matrix-valued functions F : C —>■ satisfying /jjTrace[F(/a))F(/a))^]£/a) < 0 ° 

whose components f\j are analytic in the open right half plane. For H,G € J^(iR), 
we dehne the inner-product 

(H,G)^,= trace (^{iO})G{iO})^^d(0, 

1 

with corresponding induced-norm ||H||^ = Let hnally denote by F'{X) = 

dFlds\,^X- 

The paper is organized as follows: Section [^recalls some preliminary results on 
the rational interpolation Loewner framework proposed in |[T|. Section [^presents the 
extension of these results to the single-delay case. Section derives the first order 
optimality conditions from the ^^-optimisation Problem Then, Section details 
an iterative algorithm celebrated as dTF-IRKi^ (inspired by the TF-IRKA from B). 
which allows to obtain an approximation satisfying some optimality conditions in a 
numerically efficient and memory affordable way. Finally, Section illustrates the 
proposed approach and framework on numerical examples. 

2 Realization-less interpolation 

2.1 Preliminary results in Loewner framework for rational inter¬ 
polation 

The interpolation problem, in its basic general form, is stated as follows: 

Problem 2 (General interpolation problem |[T1) Given right.- 

{(1,-, r,-, w,-) I A; e C, r,- G , w,- G ', / = 1,..., p} (4) 

and left interpolation data.- 

{(P;,1„V,)|P; G C,l; G j G C'"j = 1 , . . . , V} (5) 

construct a realization H = (E^A,B,C) of appropriate dimensions whose transfer func¬ 
tion H{s) = C{sE —A)^^B both satisfies the right.- 

//(A,)r,=w/, 1= l,...p (6) 

and the left constraints.- 

=V;, ;■= (7) 

The above problem can be solved (to obtain real-valued matrices) thanks to the 
following theorem, proposed by IT] . 

Theorem 1 (Loewner framework IQ) Given right and left interpolation data as in 
and assuming that p = V = r , the realization H = (E,A,B,C) of order r con¬ 
structed as 

£ = -L,A = -Lcr,B = y,C = VT, (8) 


^dTF-IRKA stands for delay Transfer Function Iterative Rational Krylov Algorithm. 





interpolates the right and left constraints if 

ry..^ v/r; -1,W; \i{H{li)-H{pj))rj 

ru 1 - _ Pih{H{Xi) -H{pj))rjXj 


( 9 ) 


known as the Loewner and the shifted Loewner matrices, respectively, and 

Theorem [T] allows to obtain a model H = {E,A,B,C) whose transfer function in¬ 
terpolates right and left constraints as stated in Problem]^ This has been extensively 
used for system identification from complex data obtained by a signal generator and 
for large-scale model approximation purposes llT3l . llT4l . An extension of Problem]^ 
including some derivative constraints, has also been considered to solve the model 
approximation problem 111. To this aim, the following theorem, initially stated in HI, 
provides a solution for the problem with derivatives constraints in the case where the 
right and left interpolation points are equals, i.e., Si = /r, = A, , Vi = 1,... r. 


Theorem 2 (Derivative Loevraer framework ||T||) Given a system represented by its 
transfer function H{s), r shift points {si,..., Sr} G C and r tangential directions {li,..., 1^ 
{ri,...,rr} G the r-dimensional descriptor model H = (E,A,B,C), as 

in 'll. interpolates H{s) as follows, for k = 1,... r; 




hH'{sk)rk = UH'{sk)rk, 


ifforiJ= l,...,r: 


{E)ij = 


\i{H{si)-H{sj))rj ^ 

^ J 

Si-Sj 

-liH'{si)ri i = j 


iA)ij = 


h{siH{si)-SjH{sj))rj . ^ 

^ ^ J 

Si-Sj 

-l,(s//(s))'|,.=,,.r,- i = j 


C 


[H{si)ri,...,H{sr)rr] andB = 


hH{si) 

lrH{Sr) 


( 10 ) 

( 11 ) 


In the following section the extensions of both Theorems [T] and are presented in 
the case where a time-delay realization {E,A,B,C, t) as in Q is looked for. 


3 Delay Loewner framework 

Before introducing the main result of this section, let us consider the following repre¬ 
sentation of system <11. which will be useful along the rest of the paper. 











Lemma 1 Given = (£’,A,5,C, t), its transfer function H(i{s) can be decomposed 
as: 

Hd{s)=G{f{s))e^^ (12) 

where G{s) is the transfer function of the delay-free model G = {E,A^B,C) as in 
and f{s) = se'"^. 

Proof 1 The result is straightforwardly obtained by injecting f{s) in as: 

Hd{s) = C{sE-Ae-^^)-'^B 

= clse'^E-AY'^Be'^ ( 13 ) 

= G{se^Ye^\ 

Then one extension of Theorem [T| which makes feasible the interpolation with a 
single delay descriptor system as defined in Q can be done by using f{s) as a variable 
substitution and applying the standard Loewner framework to the new transformed 
data. This first main result can be stated as follows: 

Theorem 3 (Delay Loewner framework) Let us consider p = v = r, tsM and given 
(A,,r,, W;) and the right and left interpolation data respectively, as stated in 

Assuming that f{s) = se^^ is one-to-one in the interpolation points domai^^ 
and let G = {E,A,B,C) be a realization satisfying right and left constraints from the 
data (/(Aj),ri,Wje^'^'^) and {f{pLj),lj,\je^^Y constructed with Theorem\l\ Then 
Hd = {E,A,B,C,t) satisfies the right.- 


Hd{Xi)ri = Wi, /= l,...r 

(14) 

and left constraints.- 

for the given right and left interpolation data. 

(15) 

Proof 2 The result for the right constraints (|14[) is obtained as follows: first note that 
if the delay-free model G{s) satisfies the right constraints for (/(A,),r,', then 

one obtains: 

G(/(A,))r; = w,c-^'^ (16) 

then, it equivalently follows that: 

G{f{Xi))e^‘^ri = w; 

(17) 

and by invoking Lemma^ we obtain the result: 


Hd{Xi)ri = w,-. 

(18) 


The left data constraints is similarly obtained. 

Theorem 1^ provides a method to construct a model Hd = {E,A,B,C,x) whose 
transfer function H^s) = C{sE — Ae^^''^Y^B interpolates given right and left con¬ 
straints. This is possible by noticing that the problem can be rewritten as right and 
left interpolation constraints for the delay-free for which a realization is obtained by 
the standard Loewner framework as in Theorem [T] A similar reasoning enables the 
generalization of Theorem|^as stated follows. 

"*This means that for any h{,h 2 £ {Ai,...,Ar} U {/ij ,...then f{hi) f ffii) if hi f- h 2 , where 
f{s) = 




Theorem 4 (Derivative delay Loewner framework) Let us consider a given system 
represented by its transfer function H{s), r shift points {ii,..., G C and r tangential 
directions {li,... ,1^} £ {ri,... ,rr} £ We assume that for all k m, 

fi^k) fi^m), where f{s) = se^^ (f is one-to-one in the interpolation points domain). 
The r-dimensional single delay model H = {E,A,B,C,x), as in Q. interpolates H(s) 
as follows, for k= 1,..., r; 

H{sk)rk = H{sk)rk, \kH{sk) = (19) 

\kH'{sk)bk = \kH'{sk)rk, (20) 

if only if the r-dimensional delay-free model G = {E,A,B,C) is constructed with the 
derivative Loewner framework as in Theorem^^or the shift points: 

{ai,...,ar) = {fisi),...,f{sr)), 

and the transfer function evaluation: 

(G(aiG((J,)) = {H{s,)e-^^\...,H{sr)e-^^^) 

and the derivative transfer function evaluation: 

{G'iai),...,G\ar)) = {Fi,...,Fr) 

where for i — 1 ,... r; 

G'iOi) = FiH{si),Hfsi),Si)=Fi 

Proofs First, one can note that a single delay descriptor system can be expressed as 

H^{s) = C{se-'"'E - = G{f{s))e''^ (25) 

where G{s) is a descriptor system whose representation is {E,A,B,C) and f{s) = 
se^^. Thus one can use Loewner matrices to construct the realization of system G{s) 
for the shift points ((Ji,...(7r) = (/(■Sl)) ■ ■ ■ )/('Sr)) whose transfer function data are 
(G((7i), ... ,G{<Jr)) = {H(si)e^^'^,... ,H{sr)e^^'''^). For the transfer function deriva¬ 
tive data, one can take the derivative of ( |25| l written as G{f{s)) — Flci{s)e^^'^ with 
respect to s as follows 


( 21 ) 

( 22 ) 

(23) 

(24) 


G(f{s))f'{s)=H',{s)e-^^-THMe-^^ 

and by solving the equation for G'{sk) one can obtain the result. 

This theorem allows to obtain a single delay descriptor system which interpolates 
any given transfer function F[{s). This can also be used in the case of data obtained 
through a signal generator, considering that the derivative is accessible as well. Appli¬ 
cations of this result can be found in section|6l 

Now that the delay interpolation framework has been established, one might be 
interested in obtaining a good interpolant in the sense of the J^-norm as formulated in 
Problem[^ We will now formulate mathematical conditions to select the optimal shift 
complex points i, and tangential directions r, and 1, . 



4 3^2 model reduction optimality conditions 


4.1 Preliminary results in M 2 model reduction optimality condi¬ 
tions 

The first-order optimality conditions for the delay-free H = [E,A,B,C^Q) Problemj^in 
terms of poles and residues are given in Theorem]^ 

Theorem 5 ( 0121 1 Assume that H and H have semi-simple poles and suppose that H 
is a r’^'-order finite-dimensional model with transfer function: 


His) 


k=\ ^ 


(26) 


^H,H € M 2 and H of the form (|^, is a local minimum of the M 2 delay-free approxi¬ 
mation problem, then the following interpolation equations hold: 

Hi-Xk)bk = Hi-Xk)bk, clHi-X^) = cjHi-Xk) (27) 

clH'i-X^)bk = clH'i-X^)h, (28) 

for all k = 1,... ,r where Xk kire the poles o/H and bj^ and q are its tangential direc¬ 
tions, respectively. 

As previously for the interpolation conditions, the ,^-optimality conditions are now 
extended to the single delay case. 


4.2 Results in single delay model reduction M 2 optimality condi¬ 
tions 

Proposition 1 Using the notation in Lemma^ X £ C is a pole of H^^is) = C{sE — 
Ae^^^)^^B if and only if fiX) G C is a generalized eigenvalue of the pair {A,E). 

Proof 4 X is pole of His) <;=^ w(A) is pole of Gis) w(A) G GiA,E), the 

spectrum of pencil (A,is). 

Now let recall that the Lambert W function is a multivalued (except at 0) function 
associating for each If^' complex branch, a complex number Wjt(z) such that 

2 = Wkiz)e'^kU) ^ kGZ. 

Consequently any TDS can be viewed as a system with an infinite number of poles. 

Now, as in the delay-free case, analog optimality conditions are derived in the case 
where the approximation model has the form ([T]). 

Theorem 6 (Delay -optimality conditions) Assume that H and have semi-simple 
poles and suppose that Hj is a f^-order single delay model whose transfer reads: 

r p.uT 

H,3)=CisE-Ae-^3-^B^Y. ( 29 ) 

f^^s-Uke 

TjCHjHd G M 2 and H is a local minimum of the M 2 approximation problem, then the 
following interpolation equations hold: 

Hi-Xk.p)bk = Hdi-Xk,p)bk, c^Hi-Xk^p) = clHdi-Xk^p) 


(30) 




( 31 ) 


clH{-Xkp)bk = clHd{-Xk^p)bk, 

for all p = 1,..., r and k gZ where, for given p, Xk.p are defined by: 

h.p = \wk{xap) (32) 

where Wk is the branch of the multivalued Lambert function. 

Proof 5 The proof is similar those of Theorem^as in nis using the infinite poles and 
residues decomposition of model Hd{s). 

Theorem Instates that the optimal model Hj of Problem if it exists, satisfies an 
infinite number of optimality conditions related with the Lambert W function and the 
general eigenvalues of {E,A). Nevertheless, given T € K, as Hd{s) = C{sE — B 

is parametrized by a finite number of variables, and it can be shown that it lives in a 
sub-manifold of dimension n(n„ + ny). This can be simply shown by noticing that there 
is a simple isomorphism between {E,A,B,C,t) and {E,A,B,C,Q) and the last one is 
parametrized by «(«„ +ny) variables as it can be seen in insiini. All the optimality 
conditions cannot be achieved in the general case. However, as stated in the following 
proposition, in a given particular case, the infinite optimality conditions fall into an 
equivalent finite number of relationships. 

Proposition 2 Assuming the same hypotheses ofTheoreiri^about H and H. Moreover 
if model H = {E,A,B,C, t), then the infinite optimality conditions ofTheorem^can be 
resumed to a finite number of optimality conditions as follows for p = 1,... r .• 

H{—X\,p)bk=H{—Xi,p)bk, clH{—Xyp) = clH{—Xyp) (33) 

clH{-Xi^p)bk = clH{-Xup)bk, (34) 

where X\^p = lW\{tap) andW\ is the evaluation of the Lambert function along its 
branch. 

Proof 6 One have to prove that the finite conditions (l^-@ imply This is 

possible due to the fact that: 

fiXvp) = = Xk,pe^^'’^ = fiXk,p), VkGZ 

Thus, using the decomposition given in Lemma^ it can be shown that: 

H{Xk,p) = G{f{Xk,p))e^^>^.p 

and finally 

H{Xk^p)hk = CH{Xi,p)hk = CH{Xi,p)hk = H{Xk,p)hk, 

yk G Z, where C = phg reasoning is analog for the right and derivative 

constraints, which concludes the proof. 

Now that the optimality conditions have been derived, next section is dedicated to 
the derivation of an algorithm, based on TF-IRKA Q, denoted dTF-IRKA for delay 
Transfer Function Iterative Rational Krylov Algorithm, which allows to obtain a sub- 
optimal model of the form ([T]), satisfying n(n„ +ny) optimality conditions. 


5 Delay TF-IRKA algorithm 


The algorithm proposed in this section permits to derive a system which satisfies the 
optimality conditions for r complex points. The idea behind is based on TF-IRKA 
121 which finds a model satisfying the optimality conditions in (|27li-(|20li using a fixed 
point iteration. For each iteration the new shift points will be the poles located in the 1^' 
branch of the Lambert function, only. This algorithm is celebrated as delay TF-IRKA 
(or dTF-IRKA) and is summed up as follows: 


Algorithm 1 dTF-IRKA 

1: Initialization: transfer function H{s), approximation order r S N*, 

<7** = {(jf,..., O’®} € C initial interpolation points and tangential directions 
{hi,..., br} € C""^' and (ci,...,c^j G ', W\ the first branch of the Lambert 
function. 

2: while not convergence do 

3: Build (E, a, B, C, t) using Theoremj^ 

4: Solve the generalized eigenvalue problem in yf^ and for i= 


f)*£Wxf^ = 5 ,,, 


5: 5ef, for ! = 1,... ,r 


‘ (<:+!) 

C,’ 


•(— 


e«xf) 


(35) 


(36) 


6 

7 
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end while 
Ensure conditions ([3^-(OT 
Build H = (E, A, B^C, T 


are satisfied. 


If the algorithm converges, the approximation model will satisfy optimality condi¬ 
tions given by ([33]l-([34ll and will therefore be suboptimal. The dTF-IRKA then allows 
to obtain good (in the sense of the metric given in Problem [TJ shift points and tangen¬ 
tial directions for which the interpolation problem will lead to a good approximation 
model. 

As a remark, one should note that the Lambert function evaluated in the F' branch 
can sometimes associate a real number to a complex one. In this way, the shift points 
might not be a closed set (by conjugation) and the obtained single delay interpolation 
model will not have a real representation. To avoid this, one should enforce at each 
iteration the shift points to be closed by conjugation. 

6 Applications 

This section is dedicated to the application of both methods proposed in Sectionsj^and 
1^ namely, the delay model interpolation and optimal model approximation. We 
will emphasize the potential benefit and effectiveness of the proposed approach. 






6.1 Example 1: rational interpolation 

Let us consider a dynamical model governed by the following delay model H G M 2 
whose transfer function is given as 


H{s) = 


2s+1.3e-" 


(37) 


First, model ( [J7] i (which is obviously of order 2) is approximated by a delay-free 
model of order r = 2 using the TF-IRKA (Figure [T] green dashed dotted curve). It 
is also interpolated using the delay Loewner framework with derivatives as stated in 
Theorem whose delay is set to T = 1, at the shift points si = 0.1 and S 2 = 1 (Figure 
[T] red dashed thick curve). All the results are reported on Figure and compared to 
the original model H[s) (solid blue line). 


Bode Diagram 



Figure 1: Bode diagram of original model (blue solid line), model of order r = 2 ap¬ 
proximated with TF-IRKA (green dashed dotted curve) and delay interpolation model 
using Theorem|^of order r = 2 (red dashed line). 

Figure shows that model dehned in ( |J7] i is well interpolated by a delay model 
obtained by Theorem]^ for any interpolation points. Indeed, since the transfer func¬ 
tion ( |T7| ) has a realization of the form ([T]) of order 2 it can be reconstructed using the 
Theorem]^ Figureshows quite similar results but where TF-IRKA has targeted an 
order r = 4. 

This specihc example clearly emphasizes the fact that, if the original model is a 
delay model, the counter part of obtaining a good delay-free approximation (e.g., using 
TF-IRKA) is to increase the approximation order (here the original model of order 2 
must be approximated with an order 4 to well recover the frequency behaviour). As 
illustrated on Figure]^ even with an order r = 4, the delay-free model cannot perfectly 
recover the original infinite dimensional model, while the delay model (obtained by 
Theorem]^ provides perfect matching (subject to numerical machine precision errors). 
On the other hand, the proposed delay Loewner framework allows to hnd an exact 
realization. 














Bode Diagram 



Figure 2: Bode diagram of original model (blue solid line), model of order r = 4 ap¬ 
proximated with TF-IRKA (green dashed dotted curve) and delay interpolation model 
using Theorem|^of order r = 2 (red dashed line). 


6.2 Example 2: optimal approximation and method scalability 

Let us now consider the SISO Los-Angeles Hospital model extracted from the COMPleib 
library ifT^ whose order isn — 48, denoted M^huUd = C{sl 4 ^ — S In order to 

fit the framework proposed in this paper, a delay model is constructed by injecting an 
internal delay T = 0.01 to all states, i.e., ^deiay = ^(■5^48 — This last transfer 

function is firstly interpolated on the basis of realisation of order r = 10 by applying 
the delay Loewner framework from Theorem using 10 real shift points logarithmi¬ 
cally spaced from 0.1 to 1. Then, an approximation is obtained using the dTF-IRKA 
algorithm proposed in Section]^ Figure [^compares the Bode plot of these models. 
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Figure 3; Singular value frequency response diagram of original model (blue solid 
line), approximation error with model of order r = 2 obtained with TF-IRKA (green 
dashed dotted curve) and approximation error of the delay interpolation model using 
Theorem]^ of order r = 2 (red dashed line). 



Figure 4; Bode diagram of original model (blue curve), delay Loewner interpolation 
model of order r = 10 (green dashed curve) and dTF-IRKA of order r = 10 (red dashed 
curve). 

As clearly shown on Figure 1^ the proposed dTF-IRKA allows to obtain shift 
points and tangential directions for which the interpolated delay model is much more 
accurate than the approximation using random shift points. This shows the scalability 
of the proposed approach for larger models. 






















7 Conclusions and perspectives 


In this paper, the problem of interpolating and approximating any dynamical model 
(provided its transfer function or its evaluation at given points) by a single time-delay 
finite dimensional one is analysed. Firstly, we present an extended framework which 
generalizes the Loewner one HI to the case where the interpolant is a single time-delay 
model. Then, as a second contribution, the J^-optimality conditions are derived to 
solve Problem[T] leading to an infinite set of conditions. Finally, an algorithm, denoted 
dTF-IRKA, allowing to obtain a model which satisfies a finite number of optimality 
conditions is developed and successfully applied to some numerical example; ^ 

One weakness of the proposed method, is the fact that one should know in advance 
the delay value T. Future works will investigate this issue by taking into consideration 
the delay as an optimization variable in the 3^% optimization problem. The extension 
to the multiple delay case will also be addressed in future works. 
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